
function y=prediction_MU_basic_rewt_comm(work_dir_path, beta, theta_MU, pi, R)

eval_data_source=[work_dir_path '\EVAL_DATA']; 

%data cycle
%rep=1: in-sample analysis
%rep=2:(R+1): leave-one-out bootstrap 

for rep=1:(R+1)   
    
    %Load evaluation data set
    data=load([eval_data_source num2str(rep) '.TXT']);
    
    dim=size(data);
    K=dim(2)-1; %the number of raw covariates (excluding constant, transforms, etc)
    N_obs=dim(1); %number of observations used for evaluation
    
    y=data(:,1);
    y=y-(y==0); %y=0 --> y=-1
    x=data(:, 2:K+1);
    
    c=cutoff1(x); %N_obs x 1 vector
    b=b1(x);      %N_obs x 1 vector
    
    S1 = b.*( y+1-2*c ); %N_obs x 1 vector
        
    % 1. Evaluation under perfect foresight
    %***************************************
    
    pred_PF=y;
    
    [S_PF(rep, 1) TP_PF(rep, 1) TN_PF(rep, 1) FP_PF(rep, 1) FN_PF(rep, 1)]...
        =evaluate_pred_rewt_comm(y, pred_PF, S1, pi);
    
    
    % 2. Evaluation of the logit (ML) models
    %*****************************************
    
    x_l=x_logit1(x);
    index_logit=x_l*beta(rep,:)'; %N_obs x 1
    prob_logit_at_eval_data=Lambda(index_logit); %N_obs x 1
    pred_logit=(prob_logit_at_eval_data>c)-(prob_logit_at_eval_data<=c);
    
    [S_ML(rep,1) TP_ML(rep,1) TN_ML(rep,1) FP_ML(rep,1) FN_ML(rep,1)]...
        =evaluate_pred_rewt_comm(y, pred_logit, S1, pi);
    
    %Efron and Tibshirani (1997)
    if rep==1 %in sample
        
        gammahat_ML=(  (1-pi)*sum(  b.*c     .*(pred_logit==-1) )...
            +     pi  *sum((-b.*(1-c)).*(pred_logit==-1) )  )/N_obs
        
    end
    
    % 3. Evaluation of the models fit by MU
    %***************************************
    
    x_m=x_m1(x);
    prob_MU_at_eval_data=m1( x_m, theta_MU(rep,:) ); %N_obs x 1
    pred_MU=(prob_MU_at_eval_data>c)-(prob_MU_at_eval_data<=c); %N_obs x 1
    
    [S_MU(rep,1) TP_MU(rep,1) TN_MU(rep,1) FP_MU(rep,1) FN_MU(rep,1)]...
        =evaluate_pred_rewt_comm(y, pred_MU, S1, pi);
    
    %Efron and Tibshirani (1997)
    if rep==1 %in sample
         
        gammahat_MU=(  (1-pi)*sum(  b.*c     .*(pred_MU==-1) )...
            +     pi  *sum((-b.*(1-c)).*(pred_MU==-1) )  )/N_obs
         
    end
    
    % 4. Evaluation of the WRA decision rule
    %***************************************
    
    pred_WRA=1*(sum(x,2)>=3);
    pred_WRA=pred_WRA - (pred_WRA==0); %0 --> -1
    
    [S_WRA(rep, 1) TP_WRA(rep, 1) TN_WRA(rep, 1) FP_WRA(rep, 1) FN_WRA(rep, 1)]...
        =evaluate_pred_rewt_comm(y, pred_WRA, S1, pi);
    
end %data cycle

%********************************************************************
%          AVERAGING OVER BOOTSTRAP SAMPLES (if applicable)
%********************************************************************

if R>0

    G=R;
    g=[2:R];

    %1. Perfect foresight
    av_S_PF_out  = sum( S_PF(g, 1) )/G;
    av_TP_PF_out = sum( TP_PF(g, 1) )/G;
    av_TN_PF_out = sum( TN_PF(g, 1) )/G;
    av_FP_PF_out = sum( FP_PF(g, 1) )/G;
    av_FN_PF_out = sum( FN_PF(g, 1) )/G;

    %2. Logit models
    av_S_ML_out = sum(  S_ML(g,1) )/G;
    av_TP_ML_out = sum( TP_ML(g,1) )/G;
    av_TN_ML_out = sum( TN_ML(g,1) )/G;
    av_FP_ML_out = sum( FP_ML(g,1) )/G;
    av_FN_ML_out = sum( FN_ML(g,1) )/G;

    %3. MU models
    av_S_MU_out = sum(  S_MU(g,1) )/G;
    av_TP_MU_out = sum( TP_MU(g,1) )/G;
    av_TN_MU_out = sum( TN_MU(g,1) )/G;
    av_FP_MU_out = sum( FP_MU(g,1) )/G;
    av_FN_MU_out = sum( FN_MU(g,1) )/G;

    %4. WRA decision
    av_S_WRA_out = sum( S_WRA(g, 1) )/G; 
    av_TP_WRA_out = sum( TP_WRA(g, 1) )/G;
    av_TN_WRA_out = sum( TN_WRA(g, 1) )/G;
    av_FP_WRA_out = sum( FP_WRA(g, 1) )/G;
    av_FN_WRA_out = sum( FN_WRA(g, 1) )/G;

end

%********************************************************************
%         COMPUTING .632+ BOOTSTRAP ESTIMATES (if applicable)
%********************************************************************


if R>0

%Computing weights

%Relative overfitting measure
Rhat_ML=( S_ML(1) - av_S_ML_out )/(S_ML(1)-gammahat_ML)
Rhat_MU=( S_MU(1) - av_S_MU_out )/(S_MU(1)-gammahat_MU)

%Weights on bootstrap 
w_out_ML=.632/(1-.368*Rhat_ML)
w_out_MU=.632/(1-.368*Rhat_MU)

%Logit models
S_ML_632p  = w_out_ML*av_S_ML_out + (1-w_out_ML)*S_ML(1);
TP_ML_632p = w_out_ML*av_TP_ML_out + (1-w_out_ML)*TP_ML(1);
TN_ML_632p = w_out_ML*av_TN_ML_out + (1-w_out_ML)*TN_ML(1);
FP_ML_632p = w_out_ML*av_FP_ML_out + (1-w_out_ML)*FP_ML(1);
FN_ML_632p = w_out_ML*av_FN_ML_out + (1-w_out_ML)*FN_ML(1);

%MU models
S_MU_632p  = w_out_MU*av_S_MU_out + (1-w_out_MU)*S_MU(1);
TP_MU_632p = w_out_MU*av_TP_MU_out + (1-w_out_MU)*TP_MU(1);
TN_MU_632p = w_out_MU*av_TN_MU_out + (1-w_out_MU)*TN_MU(1);
FP_MU_632p = w_out_MU*av_FP_MU_out + (1-w_out_MU)*FP_MU(1);
FN_MU_632p = w_out_MU*av_FN_MU_out + (1-w_out_MU)*FN_MU(1);

end

%**************************************************************************
%                           REPORT RESULTS
%**************************************************************************

%NOTE: The function s( ) formats the numerical output (rounds to three
%decimals)

%DISPLAY IN-SAMPLE RESULTS
%**************************

disp(' ')
disp('\begin{table}[!p]')
disp('{\footnotesize')
disp(['\captionstyle{center}\caption{$n=370$; ' '$\tau=' num2str(pi) '$}'])
disp('\begin{center}\label{tbl: ***}')
disp('\begin{tabular}{lccccccc}')
disp('\hline\hline')
disp('Method & TP & TN & FP & FN & Correct & E[U] & $\frac{E[U]}{E[U_{PF}]}$\\')
disp('\hline\hline')
disp('\multicolumn{8}{c}{IN-SAMPLE}\\')
disp('\hline')

%Display perfect foresight
disp(['PF &' s(TP_PF(1)) '&' s(TN_PF(1)) '&' s(FP_PF(1)) '&' s(FN_PF(1)) '&' s(TP_PF(1)+TN_PF(1)) '&'...
    s(S_PF(1)) '&' s(100*S_PF(1)/S_PF(1)) '\\[-5pt]'])

%Display logit model
disp(['ML &' s(TP_ML(1)) '&' s(TN_ML(1)) '&' s(FP_ML(1)) '&' s(FN_ML(1)) '&' s(TP_ML(1)+TN_ML(1)) '&'...
    s(S_ML(1)) '&' s(100*S_ML(1)/S_PF(1)) '\\[-5pt]'])

%Display MU model
disp(['MU &' s(TP_MU(1)) '&' s(TN_MU(1)) '&' s(FP_MU(1)) '&' s(FN_MU(1)) '&' s(TP_MU(1)+TN_MU(1)) '&'...
    s(S_MU(1)) '&' s(100*S_MU(1)/S_PF(1)) '\\[-5pt]'])

%Display WRA 
disp(['WRA (weed if  $\ge 3$, av.)&' s(TP_WRA(1)) '&' s(TN_WRA(1)) '&' s(FP_WRA(1)) '&' s(FN_WRA(1)) '&' s(TP_WRA(1)+TN_WRA(1)) '&'...
    s(S_WRA(1)) '&' s(100*S_WRA(1)/S_PF(1)) '\\[-5pt]'])

disp('\hline')
disp('\end{tabular}')
disp('\end{center}')
disp('\smallskip')
disp('\footnotesize}')
disp('\end{table}')
disp(' ')
disp(' ')

%DISPLAY OUT-OF-SAMPLE RESULTS (if applicable)
%*********************************************

if R>0
    
    disp(' ')
    disp('\begin{table}[!p]')
    disp('{\footnotesize')
    disp(['\captionstyle{center}\caption{Bootstrap samples=$' num2str(R) '$; ' '$n_{est}\approx 0.632\times 370$; ' '$\tau=' num2str(pi) '$}'])
    disp('\begin{center}\label{tbl: ***}')
    disp('\begin{tabular}{lccccccc}')
    disp('\hline\hline')
    disp('Method & TP & TN & FP & FN & Correct & E[U] & $\frac{E[U]}{E[U_{PF}]}$\\')
    disp('\hline\hline')
    disp('\multicolumn{8}{c}{OUT-OF-SAMPLE}\\')
    disp('\hline')

    %Display logit model
    disp(['ML &' s(TP_ML_632p) '&' s(TN_ML_632p) '&' s(FP_ML_632p) '&' s(FN_ML_632p) '&' s(TP_ML_632p+TN_ML_632p) '&'...
        s(S_ML_632p) '&' s(100*S_ML_632p/S_PF(1)) '\\[-5pt]'])

    %Display MU model
    disp(['MU &' s(TP_MU_632p) '&' s(TN_MU_632p) '&' s(FP_MU_632p) '&' s(FN_MU_632p) '&' s(TP_MU_632p+TN_MU_632p) '&'...
        s(S_MU_632p) '&' s(100*S_MU_632p/S_PF(1)) '\\[-5pt]'])

    disp('\hline')
    disp('\end{tabular}')
    disp('\end{center}')
    disp('\smallskip')
    disp('\footnotesize}')
    disp('\end{table}')
    disp(' ')
    disp(' ')

end

   
   



